function [P, n] = RengeKuttaMethodWithEps( Fxy, y0, a, b, eps )
n0 = ceil((b - a) / sqrt(eps));
P1 = RungeKuttaMethod(Fxy, y0, a, b, n0);
P0 = RungeKuttaMethod(Fxy, y0, a, b, n0 * 2);
curr_eps = ComputeDeltaForPartialFunction(P1, P0);
while (curr_eps > eps)
    n0 = n0 * 2;
    P1 = P0;
    P0 = RungeKuttaMethod(Fxy, y0, a, b, n0 * 2);
    curr_eps = ComputeDeltaForPartialFunction(P1, P0);
end;
P = P1;
n = n0;
end

